% Example,
%   % Read two greyscale images of Lena
   %Imoving=imread('images/modtest1.png'); 
   %Istatic=imread('images/modtest1.png'); 
   %Istatic=imread('images/modtest1-15d.png');
   Imoving=imread('images/lenag1large.png'); 
   Istatic=imread('images/lenag2large-15d.png');
   %Imoving=imread('images/lenag1.png'); 
   %Istatic=imread('images/lenag3.png');
   %Imoving=imread('images/lena_deform.png'); 
   %Istatic=imread('images/lena_regular.png');
   %figure, imshow(Imoving);
   %figure, imshow(Istatic);
% 
%   % Register the images
   [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,struct('Similarity','p','Registration','NonRigid'));
%    [Ireg,Bx,By,Fx,Fy] = affineregister_images(Imoving,Istatic,struct('Similarity','p','Registration','NonRigid'));
%
%   % Show the registration result
   %imwrite(Ireg,'result_DemonswithAF.png');
   %imwrite(Ireg,'result_AFbuildin2.png');
   %imwrite(Ireg,'result_standardDemonsbox.png');
   %[X,Y] = ndgrid(linspace(0,1,size(Imoving,1)),linspace(0,1,size(Imoving,2)));
   [X,Y]=ndgrid(0:(size(Imoving,1)-1),0:(size(Imoving,2)-1));
   
   figure,
   subplot(2,2,1), imshow(Imoving); title('moving image');
   subplot(2,2,2), imshow(Istatic); title('static image');
   subplot(2,2,3), imshow(Ireg); title('registerd moving image');
   
   %%
   %draw defor field
   hold on 
   contour( X+Bx, 10, 'color', 'r' );
   contour( Y+By, 10, 'color', 'r' );
   hold off
   %%
   
%   % Show also the static image transformed to the moving image
   Ireg2=movepixels(Istatic,Fx,Fy,[], 3);
   
   subplot(2,2,4), imshow(Ireg2); title('registerd static image');
%
%  % Show the transformation fields
   figure,
   subplot(2,2,1), imshow(Bx,[]); title('Backward Transf. in x direction');
   subplot(2,2,2), imshow(Fx,[]); title('Forward Transf. in x direction');
   subplot(2,2,3), imshow(By,[]); title('Backward Transf. in y direction');
   subplot(2,2,4), imshow(Fy,[]); title('Forward Transf. in y direction');
%
% % Calculate strain tensors
   E = strain(Fx,Fy);
% % Show the strain tensors
   figure,
   subplot(2,2,1), imshow(E(:,:,1,1),[]); title('Strain Tensors Exx');
   subplot(2,2,2), imshow(E(:,:,1,2),[]); title('Strain Tensors Exy');
   subplot(2,2,3), imshow(E(:,:,2,1),[]); title('Strain Tensors Eyx');
   subplot(2,2,4), imshow(E(:,:,2,2),[]); title('Strain Tensors Eyy');
   
   disp(['SAD = ' num2str(sum(sum(abs(Istatic-Ireg))))]);